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Nuclear dynamics at molecule-metal interfaces: A pseudoparticle perspective 
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We discuss nuclear dynamics at molecule-metal interfaces including non-equilibrium molecular 
junctions. Starting from the many-body states (pseudoparticle) formulation of the molecule-metal 
system in the molecular vibronic basis, we introduce gradient expansion in order to reduce the 
adiabatic nuclear dynamics (that is, nuclear dynamics on a single molecular potential surface) into its 
semi-classical form while maintaining the effect of the non-adiabatic electronic transitions between 
different molecular charge states. This yields a set of equations for the nuclear dynamics in the 
presence of these non-adiabatic transitions, which reproduce surface hopping formulation in the 
limit of small metal-molecule coupling (where broadening of the molecular energy levels can be 
disregarded) and Ehrenfest dynamics (motion on the potential of mean force) when information on 
the different charging states is traced out, which is relevant when this coupling is strong. 

PACS numbers: 


I. INTRODUCTION 

The coupled electronic-nuclear dynamics in molecules 
positioned at and interacting with metal interfaces 
presents a fundamental problem that stems from the fact 
that in such systems the usual timescale separation be¬ 
tween electron and nuclear dynamics does not necessar¬ 
ily holds. The problem has gained renewed interest in 
the context of nuclear dynamics in molecular conduc¬ 
tion junctions, a presently active field of research due 
to its fundamental and applicational importance^— Ex¬ 
perimental measurements of inelastic electron tunneling 
spectroscopy)^^— and more recently fluorescenceiip— and 
Ramaiii^^— spectroscopies serve as tools capable of pro¬ 
viding information on presence of the molecule in the 
junction and extent of heating of the device. Description 
of transport)2S heating)^ instabilities )^i2^ and current 
(and light) induced chemistry^ i^^i^^ in junctions often re¬ 
quire quantum-mechanical description beyond the Born- 
Oppenheimer approximation. 

In junctions, electron transition events between 
molecule and contacts result in coupling between dif¬ 
ferent adiabatic potential surfaces, resulting in non- 
adiabatic molecular dynamics (NAMD). NAMD plays 
important role in many chemical dynamics processes, 
ranging from surface chemistry to spectroscopy, radia¬ 
tionless electronic relaxation, photochemistry, and elec¬ 
tron transfer— Considerations of non-adiabaticity are 
particularly important fro molecules that exchange elec¬ 
trons with metal or semiconductor substrates because the 
rate of this exchange can be smaller or larger than char¬ 
acteristic nuclear timescales. Consequently, NAMD can 
drastically influences the response of molecular junctions, 
and can dominate the transport behavior associated with 
many interesting phenomena ranging from current in¬ 
duced chemistry to molecular motors—^ 

Full quantum-mechanical solution of electron- 
nuclear dynamics is possible only for relatively small 
systems— Thus one has to rely on quasi-classical 


formulations— Among them Ehrenfest dynamics^^r.^ 
and fewest switches surface hopping (SH) algorithm^i^ 
are employed most often. The latter was applied to many 
problems in the gas phase)^AL and recently also to 
molecules near metallic surfaces— From theoretical 
perspective, the Ehrenfest method can be obtained as 
an expansion around the stationary (classical) solution 
of the quantum electron-nuclei problem— >2^ Originally 
surface hopping algorithm was formulated in an ad hoc 
manner— Later work has discussed its relation to the 
quantum-classical Liouville equation—.^ Such con¬ 
siderations are not readily suitable for molecule-metal 
systems that are characterized by frequent exchange 
of electrons between a molecule and an electronic 
continuum as well as broadening of the molecular levels. 
Here we focus on this type of systems. 

In the absence of molecule-metal interaction, the 
molecule is presented in terms of its many-body molec¬ 
ular states that are usually described within the Born 
Oppenheimer (BO) approximation. The latter is based 
on the assumption that nuclei are slow relative to the 
electronic dynamics. In the other extreme limit of strong 
molecule-metal interaction, where the molecule-metal 
electron exchange is also fast relative to the nuclear dy¬ 
namics, the BO approximation is set with respect to 
hybrid molecule-metal electronic states— In the inter¬ 
mediate situation of weak but non-vanishing molecule- 
metal electron exchange coupling, the BO approximation 
breaks down and the system dynamics includes transi¬ 
tions between electronic states of different charges that 
take place on timescale of the nuclear dynamics. The en¬ 
suing dynamics can be described in the basis of the BO 
states of the isolated molecule or in the BO states of the 
strongly coupled molecule-metal complex by incorporat¬ 
ing surface hopping (SH) events into the corresponding 
nuclear dynamics. In either case, the equations of motion 
used to describe the mixed quantum-classical dynamics 
have been postulated rather than derived. In the present 
communication we offer a systematic derivation of the 
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equations of motion for such systems, and discuss their 
limiting behaviors: a surface hopping algorithm in the 
limit of weak molecule-metal coupling and Ehrenfest dy¬ 
namics on the potential of mean force in the limit where 
the electron exchange rate exceeds the characteristic nu¬ 
clear dynamics. 

Our starting point is the observation that the molecu¬ 
lar process under discussion is electron transfer into and 
out of the molecule which is most naturally described 
in the language of many-body molecular (here vibronic) 
states. This in turn requires an appropriate formula¬ 
tion of transport in the same language. The goal of 
this paper is to present such a derivation, which starts 
from the full quantum-mechanical description, and step- 
by-step derives equations suitable for implementation of 
the surface hopping algorithm for non adiabatic molec¬ 
ular dynamics at molecule-metal interfaces. The pre¬ 
sented derivation extends recent considerations^^— by 
taking into account hybridization (broadening) of molec¬ 
ular states with those of the metal(s) in a rigorous way 
and by providing expressions suitable for implementation 
of the algorithm in current carrying molecular junctions. 

The structure of the paper is as follows: after in¬ 
troducing the model in Section |n] we shortly discuss 
(Section |TTT| the pseudoparticle non-equilibrium Green’s 
function (PP-NEGF) methodology, which allows to for¬ 
mulate the molecular junction problem in the language 
of many-body states, and apply it to formulate ex¬ 
act equations-of-motion (EOMs) for the electron-nuclei 
model. Next, in Section 113 we consider the first or¬ 
der gradient expansion of these equations, which casts 
the nuclear dynamics in a classical form while maintain¬ 
ing the effect of the non-adiabatic electronic transitions 
on this dynamics. This yields a general formulation of 
the non-adiabatic dynamics at molecule-metal interfaces 
with both optical (intra-molecular) and charge transfer 
events present. The resulting semiclassical EOMs can 
be used as a basis for the surface hopping treatment of 
non-adiabatic dynamics in junctions. We specialize to 
the simple model of a resonant level coupled linearly to a 
single vibration in Section |V] in order to discuss connec¬ 
tion to previous work. Section IVTl concludes. 


II. MODEL 


We start form the usual representation of the system 
where both electron (e) and nuclear (n) dynamics is taken 
into account. The Hamiltonian is 

H{r, R) = f^{R) + Vnn{R) + Te(r) + Vee{r) + %n{r, R) 

( 1 ) 

where r and R stand for the coordinates of all electrons 
and all nuclei in the system, respectively. Tn (Te) is the 
kinetic energy of the nuclei (electrons) and (Ke) is 
the Coulomb interaction between nuclei (electrons), while 
Ven is the electron-nuclear attraction. Explicit expres¬ 


sions are 


fjR) = V As 

^ ^ 2Ma 

(2) 

a,b^l 

(3) 

< 

1 

II 

r 

(4) 

t^J — l ' ' 

(5) 

Ne „ 

Ven{r,R)- EE|^, 1 

a=l 2=1 I'i -^^il 

(6) 


Here Na and Ne represent the total numbers of atoms and 
electrons in the system, respectively. Here and below we 
have utilized atomic units, i.e. mg = = h = \. 

Our goal is to describe electronic and nuclear dynam¬ 
ics in a model junction that consists of a molecule M 
coupled to a number of metallic contacts K. The latter 
are free electron reservoirs each at its own equilibrium 
(i.e. characterized by temperature Tk and electrochem¬ 
ical potential To do so we (artificially) separate 

the whole system into molecular and contacts parts and 
assume that their electronic structure has been deter¬ 
mined. Nuclear dynamics is assumed to be confined to 
the molecular region only (and from now on we reserve 
R to represent the coordinates of the molecular atoms) 
with the contacts atoms treated as static. Coupling be¬ 
tween molecule and contacts is taken (as usual) to be 
single-particle operator (i.e. electron-electron interaction 
between electrons in M and K is disregarded). Below we 
take the index k to indicate both the band and the wave 
vector of an electron and use the second quantized repre¬ 
sentation of these states. The molecular subsystem will 
be treated in the language of vibronic states, which can 
be expanded in the basis of Born-Oppenheimer state a^^^° 

<Pev{r,R) =tpe{r,R)xUR) = ( 7 ) 


This yields the junction Hamiltonian in a mixed repre¬ 
sentation, where the molecule is described in terms of its 
vibronic states while the contacts are represented in the 
single-electron second quantized form, 

H = Hm + ^2 + Vk^ (8) 

K 


where 


eiV\,e2V2^M 

Hr = ^ Skclck 

kGK 


= E E 


X 


k^K eivi ,e2V2G.M 


( 9 ) 
( 10 ) 
+ H.C. 
( 11 ) 
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where cl (ck) creates (annihilates) an electron in level k 
of the contacts, Xe^vi,e2V2 = |ei'*^i)(e2'y2| is the molecular 
Hubbard (projection) operator, and 

H^v,,e2V2 = {eiVl\HM\e2V2) ( 12 ) 

= j dr J Hnir, R) ^e2V2{r, R) 

= X! / / dR^k{n)he^vArlri,R) 

x6i(r,i?)<l>e,„,(r,E) ( 13 ) 

are matrix elements for the molecular Hamiltonian and 
coupling to contact K. Here f dvM ■ ■ ■ integrates over 
electrons on the molecule, is number of electrons in 
the state |e2U2), ^e^v^ir/ri^R) indicates vibronic state 
|eiWi) with one electron, rj, less than in the state 1621^2), 
and Oi (r, R) is a single-electron operator, which (depend¬ 
ing on the problem) can include contributions from ([ 4 ]) 
or dH). 

III. METHOD 

Evaluating the dynamics of systems described by 
Hamiltonians of the type of Eq.® in terms of the many- 
body states of the isolated system, the nonequilibrium 
atomic limit^ can be treated within a number of tech¬ 
niques. Among them are the generalized quantum mas¬ 
ter equatiou)^^— projection operator;^ Hubbard^^— 
and pseudo particle (pp) 28 | 68- _74 nonequilibrium Green’s 
functions (NEGE) formulations, numerically exact renor¬ 
malization group approaches^^^— and quantum Monte 
Carlo methodologiesi^iSSr— The latter is usually too 
heavy to be utilized in realistic simulations. 

Here we use the PP-NEGE methodology in the low¬ 
est order (non-crossing) approximation (NCA). We note 
that generalization to higher orders is straightforwardi^ 
The PP-NEGE formulation is based on the introduction 
of second quantization in the space of the many-body 
system states 

( 14 ) 


where | 0 ) is vacuum state. The creation, pl^, and an¬ 
nihilation, Pf,y, operators satisfy the usual commutation 
relations of either Fermi or Bose operators depending on 
the number of electrons in the state \ev) This formulation 
generates an extended Hilbert space, in which the phys¬ 
ical subspace is defined by the normalization condition 

^PlvPev = l ( 15 ) 


The dynamical evolution of the system is expressed in 
terms of the pseudoparticle Green function, defined on 
the Keldysh contour as 

Ge^v^,e2V2{'rl,T2) = -i{TcPe^v^{Tl) pl^y^{T2)) (16) 

where Ty is the contour ordering operator and ri^2 are the 
contour variables. For our consideration it is convenient 
to represent this Green function (GF) in a different basis 
as follows 

Gei,e2(^i,n;l?2,r2) = (17) 

Xl\i.Rl)Ge^v^,e2V2{Tl,T2)xl\{R2) 

111 ,112 

where x%{R) is the vibrational wavefunction of the BO 
approximation ([ 71 ) in the isolated molecule. The retarded 
projection of the GF ([TT]), g2(i?i,ti;i?2,t2), gives in¬ 

formation on the many-body spectral function of the sys¬ 
tem 

24eie2(77i)li; 172 ,^ 2 ) = (18) 

* {G^ei,e2(R^^^G R'^R^.) — G'^^^^^{Ri,ti \ R2,t2)) 

where G“^_g2(i?i, tp i?2,0) = G))^ (i?2,0; l?i, G), while 

its lesser projection, G^^ e2^R^RG contains infor¬ 

mation on nonequilibrium distribution in the many-body 
states space of the molecule. These projections satisfy 
the usual Dyson equation. In particular, the following 
expressions are exact 


\ev) =ply\ 0 ) 
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* ^ - H^^^{Ri)G^^^^{Ri,ti; R 2 ,t 2 ) 

= '^ J dR y ds^S<_g(i?i,ti;i?,s)G“^(i?,s;i?2,i2) + I^l^^^{Ri,ti; R, s)Gf^^^{R, s', R2,t2) 

- G< ,e(^l. ^i; 'S)^“,e2 (^’ s; i?2, t2) - G^^ ^ (i?i, G i i?, s)S<e2 ^2, G)^ 


a 9 


9 G 9 t 2 , 

= <^61.62 <^(^1 — R2) 5{tl — ^2) 


G:,.e2 (^1, G; i?2, G) - E ( Ge,.e(^l, G; i?2, G)He"2 (^ 2 ) + (i?l, G; i?2, G) 


+ ^ y di^y G; -R, s)G;,^( i?, s; i? 2 , G) + G;^,,(i?i, G; i?, s)S;,^ {R, s; i? 2 , G) 


Here 


( 19 ) 


( 20 ) 


H^,eAR)= J dr^,,{r,R)HM{r,R)lPeAr,R)=Seue.(T,,{R)+VeiR)^ + de^e^R) + L,eAR) ( 21 ) 


where Tn{R) is defined in (I2|), 14 ( 1 ?) = VnniR) + Ee{R) is the adiabatic surface {Vnn{R) is defined in ([3|)), Ef.{R) is 
the electron eigenenergy: (Te(r) + 14e(r') + Vne{r,R))iie{r,R) = Ee{R)ipe{r, R), and 


dei ,62 {R) — 

fe,,eAR) = 



dril^eiir, R) 


9^62 {r, R) 
dRa 


d 



* 

dr (r, R) 


a^Ge2 (g r) 

dRl 


E4.e2(^) 


a—1 


d 


( 22 ) 

( 23 ) 


are the intra-molecular (not related to electron transfer between molecule and contacts) non-adiabatic couplings. Note 
that these will not couple between states of different charges. I] 4 j’ 4 ’“( 1 ? 1 jG; 1 ? 2 ,G) in (ITOl) and (EUl) are the lesser, 
retarded, and advanced projections of self-energy due to coupling to metallic contacts. Explicit expression for the 
latter within the NCA is^^ 


;(i?i,ri;i? 2 ,T 2 ) =i ^ y dR[ J di ?2 Ge'_e'(i?(, n; i? 2 , T 2 ) (24) 

^ E .6i ; 62 , 6 ' (1?1, l?i , n ; i?2 , i?2, T-2 ) - ,62 ;6i .61 (1?2,1?2, r2 ; i?'i, i?l, Ti )^ 


where c5 „/ (i?i, i?(, ri; i? 2 , 1 ? 9 , T 2 ) is the correlation 
between two electron transitions from the bath to the sys¬ 
tem; one at time ri with molecular electronic state going 
from e'l —>■ ei and nuclei changing their positions from 
R[ i?i, the other at time T 2 with molecular electronic 
state undergoing transformation from 62 —>■ 62 with nu¬ 
clei moving i ?2 —>■ i? 2 - 

.6i ;62 .6' (??1: -R'l: H ; i?2 , i?2, -^2 ) = 

E x:((??i)x:'!(??'i)tf^(?? 2 )x:^ 2 (^ 2 ) ( 25 ) 

V1,v[,V2,V2 

^ ^ ^ ^{e'^v'j^,eivi),k Ski'll:'^2^ 
k^K 


where 5 fc(ri,r 2 ) = —iiTcCkiji) c\{t 2 )) is the Green’s 
function of free electron in state k of the contact K. 


IV. GRADIENT EXPANSION 

Assuming slow nuclear dynamics we perform first or¬ 
der gradient expansion with respect to time and nuclear 
coordinates, keeping the electronic dynamics as purely 
quantum. Starting with the PP-NEGF EOMs (IT^ and 
(I20I) allows to keep information on the potential energy 
surface while going to quasi-classical description of the 
nuclear motion. Eollowing the standard procedure^^ we 
transfer to the Wigner variables, introducing slow (clas- 
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sical) coordinates and time 

i?i + i?2 ti + t 2 

- t = —— 

2 2 


( 26 ) 


so that f{Ri,ti;R2,t2) —>■ f{R,t;Rq,tq) (/ is an arbi¬ 
trary correlation function), and perform Fourier trans¬ 
form in the latter 


and fast (quantum) variables 

Rq = Rl — i?2 tq = tl — t 2 , 


( 27 ) 


f{R,t;p,E) = / dRq 


dtq 


f{R,t;Rq,tq) 

( 28 ) 


I 

Performing first order gradient expansion in Ea. dT^ leads to 




G< ..{R,t-p,E) 


* 'y '^ei,e 2 ;e 3 ,e 4 ( 7 ?, p) t; p, E) 

63 ,64 


E, 

63,64 


dR' / dp' 


dE' 

27 r 


E^ 

K 


Z,e.-.e,.e, E, E') G< t;p', E') 


( 29 ) 


where 


Gei ,62163,64 (^I p) — * *^ 64 ,63^62 ,6 


_ -A 

ai? dp ^ dR 


^/ei,63(-R) + 

' 1 a ■ 

-* 1 a 


a 1 - . , a \ 

ip -=7 

1 . 2 

^ 64 , 63 ( 1 ?)+ 2^ 

^/ei,e3 {R) dei ,63 (R) ' P 

^+ 2 '^ 64.63 (i?)^j 

^/ 64,62 (R) ~ 

■ ^ 1 a ■ 

- 1 ^ 


a 1 - 4 4 a \ 

ip -\ -=; 

[ 2 

^ 64 . 62 ( 1 ?)- 2^ 

^/e 4 ,e 2 {R) ^€4 ,62 {R) * P 

- + -d 64.62 (i?)^j 


is the Liouvillian superoperator of the free molecular evolution, and 

vZ^°l,e„eAR>R',t;p,p\E,E') =S{R-R')S{p-p') 6 {E - E')Y, J dRs j dp, J ^ 


( 30 ) 


(31) 


<^ 61,63 (cg^,e/;e 2 , 6 a ( 1 ^ 1 t;p, Ps,E Eg) Cg^g,^(i?, i?s, t; P,Ps,Es E)'jG, t; p,, E,) 

+ <^62.64 g,^ {Rs,t]p,, ^^^i) (c^,g3;g,,g4 (-R, -RsG; -p,Ps, Es-E)- c^_g^.g3_g^ (-R, i?6, t;p, -Ps, E - Eg)) (32) 

+ E ((cf 4 , 6 ; 63.64 (1?, 1?', -P, p', E'-E)- <3.g,g, (i?, i?', t'p, -p', E - E')) G“ g^ (i?, t; p, E) 

+ G^ 3 ,g(i?, t;p, E)( A3;62,64(^> -P', i? - i^') - .62:63 ,6 (1?, 1?', ” 1^))) 

is the dissipation superoperator due to coupling to contact K. X>^, 6 ^: 63 . 64 (7?, R', t',p,p', E, E') in Eq. d^Tl) is the higher 
order correction to the dissipation superoperator. Its action on the Green function Gf^^^^{R',t]p\ E') is 

Vl^ll.g,.SR^R'.t-p,p',E,E')G<^,^{R',t-p',E')= 5 {R-R') 5 {p-p') 5 {E-E') ^ J dR, J dp, J ^ (33) 

i (d64,63 ( {G<; c-> - c-<} G“ + ^ [c-> - c-<] + <562,64 (G^ {c-< - c->; G<} - ^ [c-< - c->] 


1 


E (G<{ 


aG< 


c->-c-<;G“}-^[c'-'-c 


aG“ 


.K> _ „K< _ 

at J dE 


r.^K< _ K>\ r:< 


) + ( {G' ; c'‘— c 


dG^ 




c^< - c^> 


aG< 


r 


where 

/f ■ f 1. = ElA + ElA rqq'i 
dp dR dt dE dR dp ^ ’ 


is the Poisson bracket. To shorten the notation we 
dropped the arguments in (l 33 )) keeping in mind that 
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the structure of the expression follows that of Eg. (1321) . 
This correction is responsible for renormalizations of 
G<,M t]p,E) similar to those discussed, e.g., in 
Ref. l86l In what follows we disregard this correction. 


By doing so we get in Ea. (|2^ usual structure of (en¬ 
ergy and momentum resolved flavor of) quantum master 
equation and avoid complications related to consistency 
of the gradient expansion procedure^i 


I 

Performing first order gradient expansion in EOM (l20l) leads to 


63,64 ^ K ^ 


where 


•^ei ,e 2 ;e 3 ,e 4 (-^1 — ^ei ,63 *^62,64 ( -^ 2 


*^ 62,64 |^/ei,e 3 (-R) + 
~ *^ 61,63 /e 4 ,62 (R)- 


ip - 


lA 

2 5i? 


^ 19 

W -\ -=; 

2 dK 


'^ 61 , 63 ( 1 ?) + 


lA 

2 ai? 


de.,e2{R) 2 /^ 


^/ei,e3(-R) dei,e3iR) 'P 
^/e4,e2 {^) ^€4 ,62 ('^) * P 


d 1 y r-nx d 

^+2'^e,.e3(i?)^ 

d 1 9 

- + -44.e3(i?)^ 


(35) 


(36) 


(37) 


is the free propagation superoperator^^i and 

«5^.e3;e3.e4 A t'^P. E) = t'p, E) + 

SA3,e2,eAR^f^P^E) = ^ E a/ ^ 

ee,e's 

(^de 2 ,e 4 Gl^^^,^{Rs,t;ps,E,){cA>^.^^^,^{R,Rs,t;p,-ps,E - E,) - Rs,t;-p,Ps, Es - E)) (38) 

+ A,e3Ge',e, “Ps, E - E^) - .g, RsG',-P,Ps, Es - E)^ 

is the dissipation superoperator for the retarded Green function due to coupling to contact K. sAAe 2 ;e 3 ,eiiR,t',P, E) 
in Eg. (1371) is the higher order correction to the dissipation superoperator. Its action on the Green function 
Gl^^,^{R,t-,p,E) is 


Sl^le2,eAR^t-.P.E)Gl^^,^{R,t-p,E) = \Y.j a / 


dEs 

27r 


(39) 


i.„. ( - O' {C*< - c''-^ i cq + ^ [c''> - e''<] (G- _ c'9< ; cq - ^ [c'9> - c>‘<] 


r 


The structure of Eg. (15^ follows that of (1551) . which al¬ 
lows to reproduce the omitted indices. In what follows 
we disregard this correction to the dissipation matrix. 
Eqs. (1551) and (1551) are the general final results of this 
paper. Next we turn to a specific simple example. 

V. SHIFTED HARMONIC OSCILLATOR 

In order to demonstrate relation to previous work we 
now consider a simple model of molecule represented by 
single level linearly coupled to a single harmonic oscil¬ 
lator. There are only two electronic states in this prob¬ 
lem, |0) and 11), corresponding to empty and occupied 


level, respectively. The matrix representing the molecu¬ 
lar Hamiltonian, Eg. (1211) . becomes in this case 

HA,e2 (R) = ^ ^ + Ue, (R)) (40) 

where 

d 2 

Ue{R) = — -h (5e,l(£^ + Ai?) (41) 

Here e is position of the electronic level and A charac¬ 
terizes the strength of coupling between the electron and 
molecular vibration (harmonic oscillator). 

We are interested in energy resolved joint probabilities 
to observe the oscillator at point R with momentum p 

















Then we get from Eg . (12^ 
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while the electron level is empty, Po{R,t;p, E), or occu¬ 
pied, Pi{R,t;p,E). These probabilities are defined as 

PQ{R,t;p,E) = + iG^Q{R,t-,p,E) (42) 

Pi{R,t;p,E)=-iG<{R,t;p,E) (43) 


J 


d d 


d 


dt ^dR dp 


dEi 




K 


X (^[1 - fK{Eio)]Ao{R-p,EQ)Pi{R,t;p,Ei) - fK{Eio)Ai{R;p, Ei)Po{R,f,P, Eq) 


dp J 


K 


X ( fKiEio)Ai{R;p,Ei)Po{R,f,P,Eo) - [1 - fK{Eio)]Ao{R;p, Eo)Pi{RA;p, Ei) 

I- 


(44) 


(45) 


Here Ew = E^ - Eq, /xiE) = + l]-i 

is the Fermi-Dirac thermal distribution in contact K, 
Tk = 27r^^g^|T4p(5(i? —£fe) is the electron escape rate 
to contact K (wide band approximation is assumed), and 
Ae{R;p,E) = —2ImGgg(i?;p, E) is the many-body spec¬ 
tral function of the system. Eg. (ITS)) . 

Similarly, from Eq. dSSl) we get (e=0,l) 


G:e{R;p,E) 



-1 


Ue{R)-KeiR-.P,E) 


(46) 


with 


:io{R-.P,Eo) = 


(47) 




K 


dEi 

27r 


fK{Eio)Gl,iR;p,E,) 


j:1,{R;p,e^) = 


(48) 


1 

2 


[ ^[l-fK{Ew)]Glo{R;P.Eo) 

LT 'J 


Eqs. (I44l) - (l46ll are the final results of this section. In 
spite of the simplification imparted by the gradient ex¬ 
pansion their numerical solution presents a difficult task. 
Further simplifications are achieved in two limits: 

(a) In the quasi particle approximation (T —>-0), when 
^ee -i0+, we have from m and (HSI) 

A,iR;p,E)=-2lmGl,iR;p,E) (49) 

2TrS{E - p^/2 - UeiR)) 

In this limit expressions (l44l) and (1451) reduce to those 
discussed in Ref. H 

(b) When the molecule metal coupling is strong the elec¬ 
tron exchange is fast relative to the characteristic nuclear 


dynamics. In this case individual molecular electronic 
states cannot be probed and only the some of their proba¬ 
bilities is meaningful. In this case information on the dif¬ 
ferent charging states becomes redundant, and summing 
Eqs. (l44l) and (l45l) we recover the Ehrenfest dynamics 


VI. CONCLUSIONS 


We have presented derivation of expressions for non- 
adiabatic molecular dynamics in junctions starting from 
the full quantum-mechanical problem. The derivation 
starts from the exact EOMs for the pseudo particle Green 
functions, describing junction’s response in the language 
of many-body (vibronic) states of isolated molecule. Gra¬ 
dient expansion effectively separates classical nuclear 
from quantum electron dynamics, yielding a Fokker- 
Planck equation which incorporates both optical (intra¬ 
molecular) and charge-transfer electron transitions as 
sources of non-adiabatic dynamics in junctions. The 
resulting equation can be viewed as the precursor of 
the surface-hopping algorithm. Indeed the surface hop¬ 
ping procedure described in Refs4i“— is obtained in the 
limit where level broadening is disregarded. We also 
show that tracing out information on adiabatic surfaces 
leads to Ehrenfest dynamics (motion on the potential 
of mean force). Our study extends previous considera¬ 
tion by accounting for molecular hybridization with con¬ 
tacts, and by introducing formulation capable of imple¬ 
menting surface-hopping algorithm in a current carrying 
(bias induced) molecular junction. At equilibrium and in 
the limit of weak molecule-contacts coupling our results 
reduce to those of Ref. IH- Development of numerical 
codes capable of implementing the scheme is a compli¬ 
cated technical problem that is left for future effort. 
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